R Notebook

#setwd("C:/Users/erussek/forage_jsp/analysis")
setwd("/Users/evanrussek/forage_jsp/analysis")
rm(list = ls())
library(tidyr)

Attaching package: ‘tidyr’

The following object is masked from ‘package:Matrix’:

    expand

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(zoo)

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

library(ggplot2)
library(ggpubr)

Loading required package: magrittr

Attaching package: ‘magrittr’

The following object is masked from ‘package:tidyr’:

    extract

library(RcppRoll)
library(knitr)

data <- read.csv('data/run4_data.csv')
head(data)

data <- data %>% ungroup() %>% bind_cols(s_num = group_indices(., subjectID))
data <- data %>% mutate(subj = factor(s_num))

travel_keys = data %>% ungroup() %>% select(travel_key) %>% unique()
clean_subj_data <- function(subj_data, travel_keys){
  travel_key_hard = travel_keys$travel_key[1]
  travel_key_easy = travel_keys$travel_key[2]
  
  subj_data <- data %>%select(round, phase, reward_obs, reward_true, lag, exit, start_reward, n_travel_steps,
                                    travel_key, subjectID, trial_num, s_num, correct_key) %>% group_by(trial_num) %>%
    mutate(press_num = 1:n(), round = as.factor(round), 
           rep_number = case_when(trial_num < 7 ~ 1, TRUE ~ 2)) %>%
    mutate(phase = replace(phase, phase == "Harvest", "HARVEST"))

  subj_data <- subj_data %>% group_by(trial_num) %>% 
    mutate(travel_key = replace(travel_key, travel_key == "", first(travel_key[travel_key != ""]))) %>%
    mutate(travel_key_cond = case_when(travel_key == travel_key_hard   ~ "HARD",
                                       travel_key == travel_key_easy  ~ "EASY")) %>%
    filter(!is.na(travel_key_cond))
}

# clean all the data and bind
datalist <- list()
for (i in 1:10){
  subj_data <- data %>% filter(s_num == i)
  subj_data <- clean_subj_data(subj_data, travel_keys)
  subj_data <- ungroup(subj_data)
  datalist[[i]] <- subj_data
}

Unequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vector

#cdata <- do.call(rbind,datalist)
cdata <- bind_rows(datalist)

# make a function to plot reward on each game... 
plot_subj_reward_v_press <- function(subj_data){
  
  s_num <- subj_data$s_num[1]
  
  # reward plots which show the threshold... 
  rew_plot <- ggplot(subj_data, aes(x = press_num, y = reward_obs, group = round))  +
    #geom_rect(data = t1_data, aes(xmin = press_num - 0.5, xmax = press_num + 0.5, ymin = -Inf, ymax = Inf, fill = round)) +
    geom_point(aes(color = phase)) + facet_grid(n_travel_steps ~ travel_key_cond) + theme(legend.position = "none") + ylim(0,110) + ggtitle(paste("subj: ", s_num))
  
 # plot_an <- annotate_figure(plot, top = paste('subj: ', s_num))
  plot(rew_plot)
                             
}
for (s in 1:22){
  subj_data <- cdata %>% filter(s_num == s)
  plot_subj_reward_v_press(subj_data)  
}

get_trial_exits <- function(thdata){
    # get harvest
  last_phase <- last(thdata$phase)
  
  # go through data.. if last phase was harvest, remove that round... 
  if (last_phase == "HARVEST"){
    last_round <- last(thdata$round)
    # get the last reward ops...
    last_reward_ob <- last(thdata$reward_obs[!is.na(thdata$reward_obs)])
    
    #if (last_reward_ob > 8){
    thdata <- thdata %>% filter(round != last_round)
    #}
  }
  
  # now select harvest data
  thdata <- thdata %>% filter(phase == "HARVEST", !is.na(reward_obs)) # find out why reward_true has so many .na
  
  ## get either last true reward observed
  return_tbl <- thdata %>% 
                  group_by(round) %>% 
                      summarise(s_num = first(s_num), last_reward = last(reward_obs),
                          trial_num = first(trial_num),
                          n_travel_steps = first(n_travel_steps),
                          travel_key_cond = first(travel_key_cond),
                          rep_number = first(rep_number)) %>% ungroup()
  return(return_tbl)
}

#sdata %>% group_by(trial_num) %>%
#  do(get_trial_exits(.))

exit_data <- cdata %>% group_by(s_num, trial_num) %>%
  do(get_trial_exits(.)) %>% ungroup() %>% mutate(subj = as.factor(s_num))

library(ggpubr)
exit_data <- exit_data %>% mutate(round_num = as.integer(as.character(round)))
for (s in 1:22){
  p <- ggplot(exit_data %>% filter(s_num == s), aes(x = round_num, y = last_reward)) + geom_point() + geom_line() +
  facet_grid(n_travel_steps ~ travel_key_cond) + theme_minimal()  + ylim(0,100) + ylab('last reward collected') + xlab('harvest number')
  
  #plot(p)
  #ggarrange(p)
  a <- annotate_figure(p, right = "# travel steps", top = "travel step difficulty", fig.lab = paste('subj: ' , s))
  plot(a)
}

geom_path: Each group consists of only one observation. Do you need to adjust the
group aesthetic?

# do some aggregating over exit thresholds... 
# just get the mean for each subject for each timepoint, collapse accross rep number...
mn_exit <- exit_data %>% group_by(s_num, n_travel_steps, travel_key_cond) %>%
  summarise(rep_exit_thresh = mean(last_reward), trial_num = mean(trial_num)) %>% 
  group_by(s_num, n_travel_steps, travel_key_cond) %>%
  summarise(exit_thresh = mean(rep_exit_thresh), trial_num = mean(trial_num)) %>%
  mutate(subj = as.factor(s_num)) %>% ungroup()

Unequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vector

# now just plot this for each subject -- draw lines to show easy - hard effect
ggplot(mn_exit, aes(x = travel_key_cond, y = exit_thresh)) + geom_dotplot(binaxis = 'y', aes(fill = subj), dotsize = 1.2) + geom_line(aes(group = subj, color = subj), size = 1.2) +
  facet_grid(.~n_travel_steps) + theme_minimal() + labs(y = 'Mean Last Reward Collected', x = 'Travel Sequence', subtitle = '# Travel Presses') #theme(legend.position = "none")

ggplot(mn_exit, aes(x = factor(n_travel_steps), y = exit_thresh)) +
  geom_dotplot(binaxis = 'y', aes(fill = subj), dotsize = 1.2) + 
  geom_line(aes(group = subj, color = subj), size = 1.2) +
  facet_grid(.~travel_key_cond) + 
  theme_minimal() + 
  labs(y = 'Mean Last Reward Collected', x = 'Travel Key Condition', subtitle = '# Travel Presses') 

  #theme(legend.position = "none")

# plot the means for each of these...
gmn_exit <- mn_exit %>% group_by(n_travel_steps, travel_key_cond) %>% summarise(gm_thresh = mean(exit_thresh), gsd_thresh = sd(exit_thresh)/sqrt(n()))

ggplot(gmn_exit, aes(x = factor(n_travel_steps), y = gm_thresh, color = travel_key_cond)) +
  geom_line(aes(group = travel_key_cond), size = 2) +
  geom_errorbar(aes(ymin = gm_thresh - gsd_thresh, ymax = gm_thresh+gsd_thresh), width = .1, size = 2) + 
geom_point(size = 4) + ylab('group mean exit threshold') + xlab('# travel steps') + labs(color = 'travel effort') + theme_minimal()

library(optimx)
library(lmerTest)
exit_model <- lmer(last_reward ~ n_travel_steps + travel_key_cond  + (1  + n_travel_steps + travel_key_cond |subj), data = exit_data,  control = lmerControl(optimizer ="optimx",optCtrl=list(method='nlminb')))
summary(exit_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: last_reward ~ n_travel_steps + travel_key_cond + (1 + n_travel_steps +  
    travel_key_cond | subj)
   Data: exit_data
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))

REML criterion at convergence: 6956.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1827 -0.5108  0.0002  0.5046  5.3563 

Random effects:
 Groups   Name                Variance  Std.Dev. Corr       
 subj     (Intercept)         290.72330 17.0506             
          n_travel_steps        0.06055  0.2461  -0.27      
          travel_key_condHARD  68.01256  8.2470  -0.25  0.26
 Residual                      87.82008  9.3712             
Number of obs: 928, groups:  subj, 22

Fixed effects:
                    Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)         54.95291    3.71001 20.77742  14.812 1.63e-12 ***
n_travel_steps      -0.24696    0.06318 19.11207  -3.909 0.000935 ***
travel_key_condHARD -5.06293    1.89212 17.91918  -2.676 0.015465 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) n_trv_
n_trvl_stps -0.309       
trvl_k_HARD -0.256  0.206

library(broom.mixed)

Registered S3 methods overwritten by 'broom.mixed':
  method         from 
  augment.lme    broom
  augment.merMod broom
  glance.lme     broom
  glance.merMod  broom
  glance.stanreg broom
  tidy.brmsfit   broom
  tidy.gamlss    broom
  tidy.lme       broom
  tidy.merMod    broom
  tidy.rjags     broom
  tidy.stanfit   broom
  tidy.stanreg   broom

td <- tidy(exit_model, conf.int = TRUE)
tdsome <- td %>% filter(term == 'n_travel_steps' | term == 'travel_key_condHARD')

tdsome <- tdsome %>% mutate(parameter = case_when(term == 'n_travel_steps' ~ 'travel step', term == 'travel_key_condHARD' ~ 'travel effort cond'))

ggplot(tdsome, aes(estimate,parameter,color = parameter)) + 
  geom_point(size = 2)  +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), size = 1.2) + theme(legend.position = "none") + theme_minimal() + ylab('fixed effect') + xlab('estimate (# points)') + ggtitle('linear mixed effects model - effect on exit threshold') + theme(axis.text.y = element_text(angle = 90, hjust = .5))

# reaction times... 


## a bit more sensible... 
lag_data <- cdata %>% select(s_num, rep_number, travel_key_cond, n_travel_steps, phase, trial_num, lag, correct_key, round, trial_num) %>% mutate(log_lag = log(lag)) %>% group_by(s_num, trial_num, round, phase) %>% slice(-1) %>% ungroup() %>% mutate(subj = as.factor(s_num))

ggplot(lag_data, aes(x=lag)) + geom_histogram(binwidth = 10) + xlim(50,300) + facet_wrap( ~ s_num)

#ggplot(lag_data, aes(x=log_lag))  + geom_histogram(binwidth = .01) + facet_wrap( ~ s_num)

filt_lag <- lag_data %>% group_by(s_num) %>% filter(log_lag < (mean(log_lag) + 3*sd(log_lag)) , log_lag > (mean(log_lag) - 3*sd(log_lag))) %>% ungroup()

ggplot(filt_lag, aes(x=lag)) + geom_histogram(binwidth = 10)  + facet_wrap(~s_num) + xlim(50,350)

ggplot(filt_lag, aes(x=log_lag)) + geom_histogram(binwidth = .02) + facet_wrap(~s_num)

NA
NA
NA

## with a run, plot the mean lag 
round_filt_lag <- filt_lag %>% 
  group_by(s_num, n_travel_steps, travel_key_cond, round, phase) %>%
  summarise(trial_num = first(trial_num), mean_lag = mean(lag), mean_log_lag = mean(log_lag)) %>% mutate(round_num = as.integer(as.character(round)))

# plot mean lag as a function of round number
lag_by_round <- round_filt_lag %>% group_by(s_num, round_num) %>% summarise(mean_lag = mean(mean_lag))

ggplot(lag_by_round, aes(x= round_num, y = mean_lag)) + geom_point() + geom_line() + facet_wrap(~s_num)

trial_filt_lag <- round_filt_lag %>% 
  group_by(s_num, n_travel_steps, travel_key_cond, phase) %>%
  summarise(trial_num = first(trial_num), 
            mean_lag = mean(mean_lag), 
            mean_log_lag = mean(mean_log_lag)) %>% ungroup() %>% mutate(s_num = as.factor(s_num))

ggplot(trial_filt_lag, 
       aes(x = factor(n_travel_steps), y = mean_lag, fill = s_num)) + 
  geom_point(aes(color = s_num), size = 2)+ 
  geom_line(aes(group = s_num, color = s_num), size = 1.2) +
  facet_grid(phase ~ travel_key_cond) + theme_minimal()

ggplot(trial_filt_lag, 
       aes(x = travel_key_cond, y = mean_lag, fill = s_num)) + 
  geom_point(aes(color = s_num), size = 2)+ 
  geom_line(aes(group = s_num, color = s_num), size = 1.2) +
  facet_grid(phase ~ n_travel_steps) + theme_minimal()

rr

group_lag <- trial_filt_lag %>% group_by(n_travel_steps, travel_key_cond, phase) %>% summarise(gm_lag = mean(mean_lag), gml_lag = mean(mean_log_lag)) %>% ungroup()

ggplot(group_lag, aes(x = factor(n_travel_steps), y = gm_lag, color = travel_key_cond)) + geom_line(aes(group = travel_key_cond), size = 2) + facet_wrap(~phase) + geom_point(size = 2)

rr # geom_errorbar(aes(ymin = gm_thresh - gsd_thresh, ymax = gm_thresh+gsd_thresh), width = .1, size = 2) + #geom_point(size = 4) + ylab(‘group mean exit threshold’) + xlab(‘# travel steps’) + labs(color = ‘travel effort’) + theme_minimal()

rr harvest_lag <- filt_lag %>% ungroup() %>% filter(phase == ) hl_model <- lmer(log_lag ~ trial_num + n_travel_steps + travel_key_cond + (1 + trial_num + n_travel_steps + travel_key_cond |subj), data = harvest_lag, control = lmerControl(optimizer =,optCtrl=list(method=‘nlminb’)))

rr summary(hl_model)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: log_lag ~ trial_num + n_travel_steps + travel_key_cond + (1 +      trial_num + n_travel_steps + travel_key_cond | subj)
   Data: harvest_lag
Control: lmerControl(optimizer = \optimx\, optCtrl = list(method = \nlminb\))

REML criterion at convergence: -318425.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.1768 -0.4915 -0.0665  0.3443 12.1230 

Random effects:
 Groups   Name                Variance  Std.Dev. Corr             
 subj     (Intercept)         3.952e-02 0.198786                  
          trial_num           6.360e-04 0.025220 -0.62            
          n_travel_steps      2.484e-05 0.004984  0.13 -0.48      
          travel_key_condHARD 8.152e-03 0.090291 -0.01  0.16 -0.37
 Residual                     3.686e-02 0.191985                  
Number of obs: 689562, groups:  subj, 22

Fixed effects:
                      Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)          5.1600968  0.0423935 20.9468771 121.719   <2e-16 ***
trial_num           -0.0060976  0.0053850 21.1720321  -1.132    0.270    
n_travel_steps      -0.0005989  0.0010634 21.0256448  -0.563    0.579    
travel_key_condHARD  0.0150356  0.0192651 21.0227579   0.780    0.444    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trl_nm n_trv_
trial_num   -0.620              
n_trvl_stps  0.130 -0.476       
trvl_k_HARD -0.007  0.159 -0.366
convergence code: 0
Model failed to converge with max|grad| = 0.00210383 (tol = 0.002, component 1)